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Abstract 

Let x = (xi,X 2 , • • • ,wn) be a vector of real numbers, x is said to possess an integer 
relation if there exist integers a; not all zero such that aiXi + a 2 x 2 + ■ • • + a n x n = 0. 
Beginning in 1977 several algorithms (with proofs) have been discovered to recover the 
a i given x. The most efficient of these existing integer relation algorithms (in terms of 
run time and the precision required of the input) has the drawback of being very unstable 
numerically. It often requires a numeric precision level in the thousands of digits to reliably 

recover relations in modest-sized test problems. 

We present here a new algorithm for finding integer relations, which we have named 
the “PSLQ” algorithm. It is proved in this paper that the PSLQ algorithm terminates 
with a relation in a number of iterations that is bounded by a polynomial in n. Because 
this algorithm employs a numerically stable matrix reduction procedure, it is free from 
the numerical difficulties that plague other integer relation algorithms. Furthermore, its 
stability admits an efficient implementation with lower run times on average than other 
algorithms currently in use. Finally, this stability can be used to prove that lelation bounds 
obtained from computer runs using this algorithm are numerically accurate. 
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MD 20715. Internet: helamanf (Dsuper.org. Bailey is with NASA Ames Research Center, 
Mail Stop T045-1, Moffett Field, CA 94035. Internet: dbaileyOnas .nasa.gov. 
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1. Introduction 

Let x = (xi,x 2 , • • • ,x„) be a vector of real numbers, x is said to possess an integer 
relation if there exist integers a, not all zero such that aiXi + a 2 x 2 + • • • + a n x n = 0. 
By an integer relation algorithm, we mean an algorithm that is guaranteed (provided the 
computer implementation has sufficient numeric precision) to recover the integers a, and 
to produce bounds within which no integer relations can exist. 

One application of an integer relation algorithm is to solve “subset sum problems,” 
wherein one determines what subset of a certain list of integers has a given sum. In other 
words, a subset sum problem is an integer relation problem where the relation coefficients 
a, are zero or one. This application is discussed in [11] and [5]. 

Another application of an integer relation algorithm is to determine whether or not 
a certain mathematical constant, whose value can be computed to high precision, is a 
root of a polynomial of degree n or less. This can be done for a constant a by setting 
x = (l,a, a 2 , • • • ,a n_1 ) and applying an integer relation algorithm to x. If any integer 
relation is found to hold (within the limits of the available machine precision), then the 
resulting a, are precisely the coefficients of a polynomial satisfied by a (to within machine 
precision). Even if no relation is found, a calculation with an integer relation algorithm can 
establish, for example, that the constant cannot possibly satisfy any polynomial of degree 
n or less whose coefficients are smaller than a bound established by the algorithm. This 
application is discussed in [4]. 

The problem of finding integer relations among a set of real numbers was first studied by 
Euclid, who gave an iterative algorithm which, when applied to two real numbers, either 
terminates, yielding an exact relation, or produces an infinite sequence of approximate 
relations. The generalization of this problem for n > 2 has been attempted by Euler, Jacobi, 
Poincare, Minkowski, Perron, Brun, Bernstein, among others. However, none of their 
iterative algorithms have been proven to work for n > 3, and numerous counterexamples 
have been found. 

The first integer relation algorithm with the desired properties mentioned above was 
discovered by one of the authors (Ferguson) and R. Forcade in 1977 [7]. In the intervening 
years a number of other integer relation algorithms have been discovered, including a 
non-recursive variant of the original algorithm [8], the “LLL” algorithm [12], the “HJLS” 
algorithm [10] (which is based on the LLL algorithm), and the “PSOS” [4] algorithm. 
These newer algorithms are significantly faster when implemented on a computer and 
require much less precision in the input x vector to recover the relation than the original 
algorithm. 

2. The HJLS Algorithm 

The HJLS algorithm is superior among these existing integer relation algorithms in 
several respects. For one thing, it has been proven that the number of iterations required 
for HJLS to recover a relation is bounded by a polynomial in n [10], whereas proofs of 
this property are lacking for the other algorithms. Further, the HJLS algorithm appears 
to be the most efficient of these algorithms in terms of its ability to recover the relation 
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satisfied by an input vector known only to limited precision. Finally, based on the authors’ 
experience, HJLS appears to require the lowest average computer time to recover a relation 
among previously existing algorithms. 

Unfortunately, the HJLS algorithm has one serious drawback: it is extremely un- 
stable numerically. Even for modest n and small relations, enormous numeric preci- 
sion is often required fcr the algorithm to work properly. For example, consider the 
17-long vector (1, a, a 2 , • ■ • ,a 16 ), where a = 3 1 / 4 — 2 1 / 4 . This vector has the relation 
(1,0,0,0,-3860,0,0,0,-666,0,0,0,-20,0,0,0,1). Although only 100 digits or so of a are 
required as input for the relation to be recovered using HJLS, the authors have found that 
computations must be performed with a numeric precision level of over 10,000 digits. 

One outward symptom of numerical failure in a HJLS run is that the program aborts 
due to the appearance of an extremely large entry in a matrix, one which cannot be rounded 
exactly to the nearest integer in the current working precision. Another outward symptom 
of numerical failure is that a computer run with HJLS produces a bound on the norm of 
possible relations that excludes a relation known to exist. In this second type of failure, 
comparison with runs us mg higher precision reveals that some matrix entries had become 
so corrupted by numerical error that all significance had been lost. 

In some trials, computer runs using the HJLS algorithm succeed “by accident” in 
recovering a relation using only moderate levels of numeric precision, whereas the relation 
would not be recovered at that point in the computer run if it were performed with higher 
precision. But such good fortune cannot be relied on, and it is just as likely that the 
relation will be missed a the point where it should be recovered. 

If one asks what level of precision is required for a modest-sized problem before the 
computer run using HJLS is identical to one using “infinite” precision, then the answer in 
many cases appears to he thousands of digits, based the authors’ experience. “Identical 
here means that all decisions are the same and all relation bounds are the same (to 8 digits 
or so) up to and including the point of recovery. Cases that require very high precision 
are infrequent for n < 10 but are quite common for larger n. There does not appear to be 
any a priori means of determining the required HJLS working precision level for a given 
problem. 

As one would expect, these very high levels of numeric precision require large amounts of 
memory and processing time, although amazingly enough, HJLS is still faster on average 
than other previously known integer relation algorithms. However, the most significant 
disadvantage of the HJLS numerical instability is that one can never know with certainty 
that a relation of a cerlain size has been excluded, because the bound that is obtained 
after running HJLS for a while might be completely corrupted with numerical error. One’s 
confidence in a bound result can be enhanced by increasing the numeric precision and 
verifying that the sequence of norm bounds is the same up to the point determined in the 
previous run, but one can never be certain of the result. 

The root cause of this numerical instability is not known, but it is believed to derive 
from the fact that H.JLS is based upon the Gram-Schmidt orthogonalization algorithm, 
which is known to be ni merically unstable [9]. 
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3. The PSLQ Algorithm 

Recently one of the authors (Ferguson) discovered a new polynomial time integer re- 
lation algorithm. This algorithm has been named “PSLQ,” since it is based on a paxtial 
sum of squares scheme like the PSOS algorithm, yet it can be efficiently implemented with 
a LQ (lower trapezoidal-orthogonal) matrix factorization. 

The PSLQ algorithm exhibits the favorable features of the HJLS algorithm, including 
its ability to recover relations in x vectors known to only limited precision, while completely 
avoiding its catastrophic instability. Although a straightforward implementation of PSLQ 
does not appear to be faster than HJLS on average, at least for modest-sized n, the 
numerical stability of PSLQ admits an efficient implementation using a combination of 
double precision and multiprecision arithmetic that does appear to be faster than HJLS on 
average for a wide range of problem sizes. Most importantly, one can show that by using 
a working precision level that is only slightly higher than that of the input data, bound 
results obtained from computer runs are reliable. 

Let x € R n be a nonzero n-tuple x = (xi,X 2 , ■ • • , x n _i,x„). Define the partial sums of 
squares, Sj , for x by s 2 = x l • Assume that x is a unit vector, so that si = |x| = 1. 

Define the lower trapezoidal n X (n — 1) matrix H x = ( /i tJ ) by 


ki,j 

- 0 

l<i<j<n— { 

(i) 

hij 

_ s ‘+ 1 
s% 

XiXi 

iA 

II 

^ . 
IA 
3 

1 

(2) 

fo i ,j 

_ 1 J 

1 < J < i < n 

( 3 ) 


Let P be the n x n matrix given by P = HH 1 , which has rank n — 1. By expanding this 
expression, it can be seen that P = /„ — x* • x, i.e., P = (pi,j) where p, t , = 1 — x 2 and 
Pij = —x,Xj for i ^ j. From this it follows that if x-m l = 0 then Pm 1 — m*. Let | • | denote 
the Frobenius norm, i.e., \A\ 2 = Yl a ]j- Then it can be seen that |P| = \H\ = \fn — 1. 

Given an arbitrary lower trapezoidal n x (n— 1) matrix H = (hij), define an associated 
n x n lower triangular matrix D £ GL(n, Z) as follows. The entries of the n x n matrix 
D = (d,j) are given from i back to d,j by means of the recursions 


dij = 0 


1 < i < j < n 

(4) 

dij = 1 

/ 

\ 

IA 

II 

IA 

5 

( 5 ) 

d^ = nint —7— Yl 


1 < j < i < n 

(6) 


\ 1 <j<k<i ) 

The function nint denotes the nearest integer function. 

The inverse of D, namely E = D -1 , can be defined with recursions similar to the above 
for E = (e,j) by 


e t,j = 0 

e i,j = 1 

® i,j ~ y ] Ci,kdk,j 

l<j<k<i 


l < i < j < n 
l < i = j < n 
1 < j < i <n 


( 7 ) 

( 8 ) 
( 9 ) 
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for e ^_ 1 back to e,^. 

Given an arbitrary louver trapezoidal n x (n — 1) matrix H = (hij) and a fixed integer 
j, 1 < j < n - 1, define the (« - 1) x (n - 1) orthogonal matrix Gj G 0(« - 1) as follows. 
If j = n - 1 set G n -i = l n - i, the m atrix id entity. Otherwise, for j < n - 1, set h h3 = a, 


hj+ i,j+i — G 


and d = \/b 2 + c 2 . The entries of the (n — 

1 ) x (n — 1 ) matrix 

Gj = ( 9i,k ) are given by 



01,* 

= 1 

1 < i < j or j + 1 < i < n 

( 10 ) 

0.7, J 

= 6 /d 


(ID 

0A7 + 1 

= —c/d 


( 12 ) 

0 j+ 1 J 

= cjd 


(13) 

0 j+ 1 J +1 

= b/d 


(14) 

0 *,fc 

= 0 

otherwise 

(15) 


Given a fixed integer j, 1 < j < n — 1 define the n x n permutation matrix to be that 
matrix formed by exchar ging the j-th and j + 1-st rows of the n X n identity matrix 7 n . 


5. One Iteration of PSLQ 

Select a constant 7 > 2/\/3 = 1.1547 ■ • Suppose we are given three matrices, H,A,B , 
where H is a n X (n — 1) lower trapezoidal matrix and A and B are n X n integral matrices, 
with B = A -1 . An iteration of the algorithm PSLQ is defined by the following three steps. 

1 . Replace H by DH . 

2. Select an integer j, 1 < j < n — 1 such that ^ f° r a ^*, 1 < f < n — 1. 

3. Replace H by RjHGj, A by B.jDA and B by BERj. 

Theorem. Fix 7 > 2/^/3 and set 6 2 = 3/4 - I/ 7 2 . Suppose some integral linear com- 
bination of the entries of x G R n is zero, so that x has an integer relation. Let M > 1 
be the least norm of an> such relation m. Normalize x so that |x| = 1 and iterate PSLQ 
beginning with the following set of three matrices: H = H x , A = B = I n . Then 

1 . A relation for x will appear as a column of B after fewer than | ?n 2 (n + l)log(Mn 2 ) 
iterations of PSLQ. 

2. The norm of such a relation for x appearing as a column of B is no greater than 

y/n\H\\BP\M. 

3. If after a number of iterations of PSLQ no relation has yet appeared in a column of 
B , then there are no relations of norm less than the bound l/\H\. 
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6. Proof of the PSLQ Algorithm 

Suppose m £ Z n is a relation for x, so that xm t = 0 with m/0, For any integral 
invertible C £ (7T(n, Z) and orthogonal Q £ 0(n — 1), 

1 < |CW| = |CPrn*| < \CP\\m\ = |CP r ||m| = |Ctf*Q||m| (16) 

Note that after some fixed iteration of PSLQ, H = AH X G, where (7 is the product of the 
Gj. Perform Step 1. Note that 0 ^ hi^ if no relation has been found. Then for any relation 
m, including one of norm M, 

1 < |m||i/| < \rn\ I £ (n-t + IJA?,- (17) 

V !<*<«-! 

Observe that |/i M | < 1 from the beginning so that 

£ ( n ~ l + [ ) h l> < £ (n - i + l)\h iti \ < £ ( n ~ i + l)l^«.«| 1/n ( 18 ) 

l<i<n— 1 1 < i < n — 1 l<*<n— 1 

Now select the integer j, 1 < j < n — 1 as in Step 2 and assume without loss of 
generality that a = \hjj\,b = |Aj+ij|,c = |Aj+ij + i|. Then 0 < b < a/2 and 0 < c < a/7. 

Set t = \/6 2 -f c 2 /a so that t < \/l — 6 2 = iy/l/4 + I/7 2 < 1 with 7 > 2/\/3- Perform Step 
3 and denote the diagonal entries of RjHGj by h;;. Then 

£ (rc-; + i)M 1/n (19) 


< 


< 


-(n - j + l)a 1/n (l - t 1/n ) + (n - j)c 1/n (r 1/n - 1) 

+ £ (n-;+i)|M 1/n 

1 <i<n— 1 


— a 1 ^ n (l — t 1 ^”) < 1 + (n — j) 
+ Y (n-i + l)\hi,i\ 1/n 

l<t<n— 1 



£ (n — *+l)|fti,| 1/n - 

l<i<n— 1 


(l- 


6 2 

27 n 2 (n + 1) 


) 2 £ (» 


a 1/n (l -* 1/n ) 


- i + l)|Ai, I | 1/i 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 


Inequality (24) follows from the fact that the expression in braces in (22) is always at least 
one. Inequality (25) follows from the three approximations 

1 - t 1/n > 6 2 /{2n) for n > 1 
i 1 a > \hi,i\ 

hn(« + l)fl 1/fl > £ (n-i + l)|fe w | 1/B 

1 l<i<n-l 
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(26) 

(27) 

(28) 



Suppose that k iterations of PSLQ have been done and that no relation has appeared. 
Then from (17) and (25) 


1 < |m| (1 




27» 2 (n + 1 ) 


) fc n 2 


(29) 


By taking logarithms, this becomes 

k < ‘-^n 2 (n + l)log |m|n 2 (30) 

o 1 

Consider the intersection of the hyperplane perpendicular to x and the convex hull of 
the 2n vertices given by i the rows of the transpose of B. Since | det B\ = 1 there are no 
relations in the interior o:' this intersection. The radius of the largest n — 1 dimensional ball 
in this intersection is the minimum of the 2" quantities i/\6H\ where 6 = (±1,- -,±1). 
This radius is bounded below by l/(y/n\H\). Consequently, as \H\ becomes smaller this 
relation-excluding ball becomes larger. The nearest candidates for relations correspond to 
columns of B. This concludes the proof of part 1 of the theorem. Part 2 of the theorem 
follows from the fact that the radius above is less than M , the norm of a smallest relation. 
Part 3 of the theorem follows from the general inequality l/\CH x Q\ < |m| with which the 
proof began, where C — A and Q = G in the case of the PSLQ algorithm. 

Note that the numerical expression 2'y/S 2 has a minimum of 8 when 7 = 2. Thus 
it would appear from this proof that 7 = 2 is the best choice. However, in practice we 
have found that smaller values of 7, while requiring more iterations, are more effective in 
recovering relations when the input vector is known to only limited precision. 


7. Computer Implementation 

The basic PSLQ algorithm can be implemented easily using ordinary floating point 
arithmetic on a computer. Using double precision (i.e., 64-bit) arithmetic, relations of 
modest height can be recovered for n up to about 6. Beyond this level, precision is quickly 
exhausted. Thus a serious implementation of PSLQ or any other integer relation algorithm 
must employ some form of multiprecision arithmetic. The authors’ implementation of 
PSLQ employs an automatic multiprecision translator [2] and a multiprecision arithmetic 
package [3], both of which were developed by one of the authors (Bailey). 

A straightforward implementation of PSLQ as described in section 5 is adequate for 
recovering relations, bu~ its run time is not competitive with HJLS, at least for small 
n. Fortunately, the stability of the PSLQ algorithm admits a multi-level implementation, 
i.e., one that utilizes two or even three levels of working precision, which results in greatly 
reduced run times. The wo-level scheme can be described as follows. To initialize, perform 
the initialization step of the standard PSLQ algorithm (i.e., normalize x. compute H x , and 
set A and B to the identity matrix) using full precision. Also set y 0 = x. Then perform 
the following three steps . 

First, convert the multiprecision arrays t/o and H as accurately as possible to the dou- 
ble precision arrays 1/0 and 11 . scaling the results to the maximum entries in //o and H , 
respectively. Set the double precision arrays A and B to the identity matrix. 
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Second, perform the basic PSLQ algorithm as described in section 5 using double 
precision arithmetic for as many iterations as possible until precision is exhausted. The 
authors have found that an adequate condition to detect exhaustion of double precision is 
that either an entry of the A matrix exceeds 10 8 or an entry of y = Byo is smaller than 
1(T 14 . 

Third, perform an “accurate update”. This is done as follows. First, replace the 
multiprecision A by A A and B by BB. Then compute K = AH X and perform a lower 
trapezoidal-orthogonal (LQ) factorization on K. This LQ factorization can be done by 
employing a multiprecision version of the LINPACK routine DQRDC [6], if desired. Note 
that when K is (uniquely) factored as LQ, where L is lower trapezoidal and Q is orthogonal, 
then since H = AH X G is lower trapezoidal, it follows that L = H and Q = G -1 . Thus 
H may be set to the lower trapezoidal result of the LQ factorization, and the orthogonal 
result of the LQ factorization may be discarded. 

Finally, set the multiprecision y = xB. If any entry of this updated y vector is zero (to 
within the multiprecision epsilon), then a relation has been detected and is contained in 
the corresponding column of the B matrix. The bound on possible relations may also be 
checked at this point using the updated H matrix. The algorithm may be continued by 
repeating the above procedure beginning with the conversion step, where t/ 0 now is set to 
the current updated y. 

One detail has been omitted in the above procedure. In some iterations, often including 
the very first iteration, the entries of the y 0 vector or the H matrix have such a large 
dynamic range that they cannot be meaningfully converted to double precision values. In 
these cases it is necessary to perform the basic algorithm using multiprecision arithmetic 
for a number of iterations until these large dynamic ranges are eliminated. 

With this implementation scheme, for n < 20 the run time is dominated by the time 
required for the double precision iterations and thus is very efficient. For larger n a three- 
level implementation can be used employing double precision, an intermediate level of 
perhaps 100 digits, and full multiprecision. With this scheme, the run time is dominated 
by the cost of double precision iterations for as large a value of n as is practical to perform 
on current scientific computer systems. 

One major advantage of the PSLQ algorithm is that it permits a very simple analysis 
affirming the accuracy of computed bound results. This can be seen from the above, 
since the H matrix can be updated directly from the original H x matrix. The formation 
of the H x matrix requires at most two operations per entry. The matrix multiplication 
H — AH X requires at most 2 n operations per entry, and the entries of A are integers. The 
LQ factorization of H using Householder transformations is known to be very stable [9], 
requiring at most 2 n operations per matrix entry. Finally, the computation of the bound 
l/|/f| is very stable, involving only the sum of squares and a division. 

In short, the H matrix, from which relation bounds are computed, is not the result of 
a long iterative computation, but instead may be computed directly from the original H x 
matrix at any point in the calculation with fewer than 4n operations per element. Thus it 
follows that the accuracy of the bound can be guaranteed provided that these operations 
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are performed with a level of precision only a few digits higher than that requiied to 
accurately represent the integer A matrix that has been developed at a given point in the 
calculation. For example, suppose that n = 25, the entries of the A matrix do not exceed 
lO 100 , and the result of each multiprecision arithmetic operation is known to be correct 
to within one bit in the last place. Then the norm bound result is guaranteed correct to 
at least eight decimal places, provided that the initial H x and updated H matrices are 
computed with at least L40 digit precision. In practice, the authors have found that a 
working precision level of only 20 digits or so beyond the size of the largest A matrix entry 
is sufficient to produce reliable bounds. 

Along this line, it should be mentioned that the inequality (16) above is true for any 
matrix norm. For example, the norm given by 


\H\ 


max 

l<t<n 


i 

l<j<n-l 


hh 


(31) 


gives better bounds than the Frobenius norm. It should also be emphasized that Part 2 
of the Theorem gives a measure of how well PSLQ is doing in regard to constructing the 
shortest integer relation. In practice, the PSLQ algorithm almost always finds the shortest 
relation. 

It can be seen from the above that the two- and three-level implementation schemes rely 
on the fact that the dynamic ranges of the y vector and the H matrix are of modest size 
for most iterations. It also appears that in most computer runs, after the first 1000 or so 
iterations, hundreds of iterations at a time may be performed in double precision between 
accurate updates. The reason for this fortunate behavior is not completely understood, but 
it may be related to a higher dimensional analog of the Kuzmin probability distribution 
of small partial quotients in continued fractions. We conjecture that the dynamic range of 
the diagonal entries of //, i.e., max; /i,., )/ min, h j:i \ , is bounded by 7 n n k after the initial 
cn 2 iterations for some fixed k and c. 


8. Performance Results 

Tables 1, 2, and 3 give performance results for the PSLQ algorithm. Both PSLQ and 
HJLS were implemented by one of the authors (Bailey) and perform high precision arith- 
metic using this author’s MPFUN multiprecision package [3]. H.ILS was “hand-coded” with 
direct calls to the MPFLN arithmetic routines. PSLQ was written in ordinary Fortran-77, 
using some Fortran subroutines from the LINPACK library [6] and a matrix multiply rou- 
tine from the LAPACK library [1], and was converted to multiprecision using this author’s 
automatic multiprecision translator [2]. 

The tests in Table 1 were constructed using pseudorandom number generators, and 
both algorithms were given the same set of ten test problems for each value of n. The 
tests were run on a single processor of a Silicon Graphics 380 workstation, which has a 
LINPACK 100 rating of 5 MFLOPS (double precision). The column headed “Dim. n” 
gives the dimension 11 of the relation vector, h is the height of the constructed relation 
coefficients: each a, except a n is chosen at random between — h and h. The column headed 
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d gives the number of significant digits in the randomly generated x vector. The columns 
headed “Min. Prec.” and “Max. Prec.” give, in digits, the minimum and maximum 
working precision levels used for the runs in a given set. The column headed “Succ.” gives 
the number of successful trials in the ten trials of a single set, i.e., the number of trials 
where the original constructed relation was recovered. The columns headed “Ave. Bound” 
give the average of the final norm bounds produced. The column headed “Ave. Time” 
gives the average CPU run time in seconds. 

The purpose of the trials in Table 1 is to compare the effectiveness of the algorithms 
in recovering relations for input vectors known to only limited precision. Thus these cases 
were chosen with values of d just sufficient so that the constructed relation is the relation 
of minimum norm. The authors found that if d is set even 5 or 10 digits below the level 
listed in Table 1, extraneous relations are recovered by PSLQ with smaller norms than 
the original constructed relation. This finding underscores the fact that PSLQ is very 
economical in the input precision required to recover a relation and is, in fact, quite close 
to the information theoretical minimum in this regard. 

The two-level variant of the PSLQ algorithm described in section 7 was employed for 
the runs in Table 1, with 7 = 2/\/3‘ 

For the HJLS runs in Table 1, the working precision level was initially set at 230 digits, 
which is somewhat greater than that for the PSLQ runs. In some of the HJLS runs, the 
calculation aborted because multiprecision numbers were encountered in the calculation 
that were larger than could be precisely rounded to the nearest integer in the current 
precision. In other runs, the HJLS algorithm terminated with the erroneous assertion that 
the norm bound had excluded the norm of the constructed relation. Whenever a HJLS 
run failed to recover the original constructed relation for one of these reasons, the working 
precision level was doubled and the run repeated, up to a precision level of 7400 digits (32 
times the original level) if necessary. 

In some of the HJLS runs, relations were recovered “by accident” — the relation would 
not have been recovered at that point in the run if the calculation had been performed 
at higher precision. In these cases the run time was taken to be the time for the run 
at the lower precision. In a number of other trials, HJLS recovered a relation, but this 
relation was an extraneous relation with a larger norm than the constructed relation. In 
two trials, numerical failure occurred even at 7100 digit precision, and these trials were 
deemed failures. Average run time and bound statistics for HJLS runs were computed only 
over runs where some relation was recovered, even though this results in lower average run 
times for HJLS than would be the case if the failures were counted. 

For those HJLS runs that required a precision level of 1850 digits or more, a version 
of the author’s HJLS program was employed that calls the advanced multiplication and 
division routines of the MPFUN package. The advanced multiplication routine employs a 
fast Fourier transform, and the advanced division routines use a Newton iteration. Usage 
of these routines resulted in much lower run times for those HJLS runs that required very 
high precision. 

The results in Table 1 show that PSLQ appears to be even more effective than HJLS 
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Dim. 

n 

h 

d 

Min. 

Prec. 

Max. 

Prec. 

Succ. 

Ave. 

Bound 

Ave. 

Time 

PSLQ 







7.971 x 10 n 



10 

12 

125 

140 

140 

10 

13 


15 

8 

130 

145 

145 

10 

9.412 x 10 7 

48 


20 

6 

150 

165 

165 

10 

1.029 x 10 6 

157 


25 

5 

160 

175 

175 

10 

7.471 x 10 4 

462 


30 

4 

200 

215 

215 

10 

6.371 x 10 3 

941 

HJLS 










10 

12 

125 

230 

460 

10 

1.265 x 10 12 

38 


15 

8 

130 

920 

7400 

9 

1.585 x 10 8 

2382 


20 

6 

150 

920 

7400 

9 

1.414 x 10 6 

2979 


25 

5 

160 

920 

3700 

4 

1.164 x 10 5 

641 


30 

4 

200 

920 

3700 

9 

1.142 x 10 4 

1126 


Table 1 : Random Relation Tests 


in recovering relations from input vectors known to only limited precision. In 9 out of a 
total of 50 trials, HJLS f ailed to recover a relation that PSLQ succeeded in recovering. In 
those cases where it did jecover a relation, HJLS is slower on average than PSLQ. In some 
individual trials, HJLS actually recovered relations faster than PSLQ, but in other trials, 
where very high numerical precision had to be employed for HJLS runs, HJLS was many 
times slower than PSLQ. The PSLQ run times varied no more than 25% within a given 
set. 

Some other results comparing PSLQ and HJLS are shown in Table 2. In these trials, the 
problem is to recover the polynomial of degree n satisfied by a , where a = 3 1 / r — 2 1 ^ s . In 
other words, the input vector x = ( 1, a, a 2 , • • ■ , a n_1 ), where n - rs + 1. Here the working 
precision was set to the level required by PSLQ or HJLS, respectively, as determined by 
some preliminary test runs, and the input vector was computed with ample precision to 
recover the relation. In this way, these tests simply compare run times, not effectiveness 
in recovering relations with limited input precision. 

The results in Table 2 demonstrate even more dramatically the erratic precision re- 
quirements and correspondingly erratic run times of HJLS. Whereas the first two problems 
required less than 100 digits working precision, and ran quite fast as a result, the last 
three problems required thousands of digits of precision, and the run times were hundreds 
of times greater than tie corresponding PSLQ run times. In fact, it appears from this 
and other runs made by the authors that the HJLS algorithm is generally unusable for 
recovering lion-trivial algebraic relations when n is greater than 25. 

A final set of results is given in Table 3. These results do not compare PSLQ with HJLS, 
but instead merely explore how large a relation can be recovered in a reasonable amount 
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of time using the PSLQ algorithm. These runs were made using the three-level variant of 
the PSLQ algorithm, as described in section 7, with 7 = 2/\/3. The trial problems are, as 
before, algebraic relations with a = 3 1//r — 2 1 ^ 9 . These runs were made on a single processor 
of a Silicon Graphics 380 system. 

The results in Table 3 show that it is feasible at present, using a RISC workstation, 
to recover relations with the PSLQ algorithm for n greater than 80 and, simultaneously, 
for relation coefficients larger than 10 13 . These figures are significantly larger than any 
with which the authors are familiar for other integer relation algorithms (compare with 
the results in [11], for example). It is worth noting that if a simple-minded, exhaustive 
search procedure had been employed to recover the polynomial satisfied by a — 3 1 / 9 — 2 1 / 9 , 
assuming only its numerical value and the fact that it is an algebraic number of degree 
not exceeding 81, with coefficients bounded by 10 14 , then it would have been necessary to 
examine more than 10 llo ° polynomials. 
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Dim. 



Working 

Max. 

Relation 

Run 


n 

T 

3 

Prec. 

Bound 

Norm 

Time 

PSLQ 

■ 

3 

3 

30 

7.414 x 10 1 

1.650 x 10 2 

1 


mm 

2 

5 

40 

3.773 x 10 2 

5.898 x 10 2 

3 


1 

3 

4 

55 

2.442 x 10 2 

6.446 x 10 2 

6 


E3 

2 

7 

70 

3.027 x 10 3 

8.248 x 10 3 

17 


1b 

3 

5 

75 

9.277 x 10 2 

2.699 x 10 3 

18 


■3 

4 

4 

75 

1.255 x 10 3 

3.917 x 10 3 

22 

HJLS 


3 

3 


7.441 x 10 1 

1.650 x 10 2 

2 


n 

2 

5 


1.835 x 10 2 

5.909 x 10 2 

3 


13 

3 

4 

■ 

1.787 x 10 2 

6.446 x 10 2 

131 


15 

2 

7 

1 

4.165 x 10 3 

8.248 x 10 3 

10480 


16 

3 

5 

H 

1.618 x 10 3 

2.699 x 10 3 

11112 


17 

4 

4 


1.390 x 10 3 

3.917 x 10 3 

11490 


Table 2: Algebraic Relation Tests 


Dim. 

n 

T 

s 

Working 

Prec. 

Max. 

Bound 

Relation 

Norm 

Run 

Time 

21 

4 

5 

110 

5.770 x 10 3 

1.173 x 10 4 

71 

26 

5 

5 

160 

3.092 x 10 4 

1.169 x 10 5 

218 

31 

5 

6 

230 

2.380 x 10 5 

6.686 x 10 5 

656 

37 

6 

6 

300 

1.034 x 10 6 

5.344 x 10 6 

1680 

43 

6 

7 

400 

1.756 x 10 7 

6.022 x 10 7 

3981 

50 

7 

7 

500 

5.871 x 10 7 

3.308 x 10 s 

8493 

57 

7 

8 

700 

2.142 x 10 9 

9.531 x 10 9 

20730 

65 

8 

8 

850 

2.028 x lO 10 

1.205 x 10 11 

46897 

73 

8 

9 

1050 

4.677 x 10 11 

3.120 x 10 12 

93368 

82 

9 

9 

1300 

5.026 x 10 12 

7.973 x 10 13 

185787 


Taole 3: Large Algebraic Relations with PSLQ 
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